arXiv: 1508.06492v2 [q-fin.CP] 7 Oct 2015 


Ninomiya-Victoir scheme: strong convergence, antithetic version 
and application to mnltilevel estimators 

A. A1 Gerbi, B. Jourdain* and E. Clement^ 


October 8, 2015 


In this paper, we are interested in the strong convergence properties of the Ninomiya-Victoir 
scheme which is known to exhibit weak convergence with order 2. We prove strong convergence 
with order 1/2. This study is aimed at analysing the use of this scheme either at each level or 
only at the finest level of a multilevel Monte Carlo estimator: indeed, the variance of a multilevel 
Monte Carlo estimator is related to the strong error between the two schemes used on the coarse 
and fine grids at each level. Recently, Giles and Szpruch proposed in [6] a scheme permitting 
to construct a multilevel Monte Carlo estimator achieving the optimal complexity O (e“^) for 
the precision e. In the same spirit, we propose a modified Ninomiya-Victoir scheme, which may 
be strongly coupled with order 1 to the Giles-Szpruch scheme at the finest level of a multilevel 
Monte Carlo estimator. Numerical experiments show that this choice improves the efficiency, 
since the order 2 of weak convergence of the Ninomiya-Victoir scheme permits to reduce the 
number of discretization levels. 


1 Introduction 


This paper is dedicated to the computation of V = E [/ (V^)], where / : M” —)■ M is a payoff 
function and Xt is the solution, at time T S to a multi-dimensional stochastic differential 
equation of the form 


dXt = b{Xt)dt + E ngApdW/, t G [0,r] 
i=i 

Xo = X. 


( 1 . 1 ) 


Here, x G M”' is the initial condition, W = (W ^,..., is a d—dimensional standard Brownian 
motion, h : M"' —M”' is the drift coefficient and : M”' —> K”, j G {1,..., d}, are the diffusion 
coefficients. 


The standard Monte Carlo method consists in estimating E [/ (Xt)] by discretizing the stochastic 
differential equation with A G N* steps and approximating the expectation using M G N* 
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independent path simulations. To be clear, the crude Monte Carlo estimator is given by 


Ycmc 


1 

M 


M 


where are independent copies of a numerical scheme with time step T/N. Under 

some regularity assumptions on the coefficients of the SDE and for a smooth payoff, it is well 
known that to ensure a root mean-square-error e, the computational cost of this method is 
where a is the order of weak convergence of the numerical scheme (see theo¬ 
rem 1 in [3]). In [ 8 ], Ninomiya and Victoir proposed a numerical scheme, achieving a = 2, 
which reduces the computational complexity compared to the Euler scheme, for which a = 1. 
In the optimal complexity O , the term 1/a is due to the bias E [/ (X-r)] —E [/ (X/Y)]. 


To remove this term, Giles introduced in [5] a multilevel Monte Carlo estimator permitting 
telescopic cancellation of the bias. The multilevel Monte Carlo estimator is built as follows 



where L G N* is the last and finest level of discretization with time-step T/ 2 ^, G 

-g vector of sample sizes at each level. Moreover, for all I G {1,... ,L}, the two 
numerical schemes X^ and X^ are simulated with the same Brownian motion. For each 
discretization level I G {0,L}, Mi independent and identically distributed path simulations 
independent from the other levels are used. The optimal complexity of this method is driven by 
the order [3 of convergence to zero of the variance V I / I Xj. ’ ) — / (Xrp ’ ) ) j which is related 
to the strong convergence order 7 of the scheme. For a Lipschitz payoff /, using the strong 
convergence properties of the scheme in the estimation of the variance, one gets /3 > 27 . For 
/? > 1, the optimal complexity is O This complexity is the same as in a simple Monte 

Carlo method with independent and identically distributed unbiased random variables. The 
condition /3 > 1 is satisfied by the Milstein scheme for which 7 = 1 . Unfortunately, to simulate 
the Milstein scheme, one needs, in general, to simulate Levy areas for which there is no known 
efficient method when the dimension of the Brownian motion d is larger than 2. Unless the 
diffusion coefficients G {I,...,(i}, are constant, the strong order of the Euler scheme is 
7 = 1 / 2 , which leads to /? = 1 and to the optimal complexity O (e~^ (log ( 7 ))^) • 


Recently, two approaches have been developed to improve the case 7 = 1/2. In [ 6 ], Giles and 
Szpruch introduced a modified Milstein scheme, with the Levy areas set to zero, and its antithetic 
version based on the swapping of each successive pair of Brownian increments in the scheme. 
Regarding the multilevel Monte Carlo estimator, at each discretization level I G {1,..., L}, on 
the finest grid, instead of using a simple scheme, Giles and Szpruch employed the arithmetic 
average of the scheme and its antithetic version as follows 



where denotes the antithetic version of the Giles-Szpruch scheme with time-step T/2^ Un¬ 
der some regularity assumptions on the coefficients of the SDE and for a smooth payoff, Giles 
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and Szpruch showed that despites 7 is equal to 1/2, f3 is equal to 2 which leads to an optimal 
complexity O • The second approach called multilevel Richardson-Romberg method and 
investigated by Lemaire and Pages in [7], fully takes advantage of the existence of a weak er¬ 
ror expansion while keeping the multilevel Monte Carlo estimator properties. The multilevel 
Richardson-Romberg estimator is a weighted version of the multilevel Monte Carlo method 
which integrates the multi-step Richardson-Romberg extrapolation developed by Pages in [10]. 
Lemaire and Pages obtained an optimal complexity O (e“^ log (^)) when (3 = 1 which improves 
the standard multilevel Monte Carlo method. When (3 > 1, the optimal complexity O is 
preserved. 


In this paper, we propose to use the Ninomiya-Victoir scheme, which is known to exhibit weak 
convergence with order 2 , on the finest grid at the last level L of a multilevel Monte Carlo 
estimator. This idea is inspired by Debrabant and Rossler [2] who suggest to use a scheme with 
high order of weak convergence on the finest grid at the finest level L of the multilevel Monte 
Carlo method. By this way, Debrabant and Rossler reduce the constant in the computational 
complexity by decreasing the number of discretization levels. 

In section 2, to derive the strong convergence order of the Ninomiya-Victoir scheme, we provide 
a suitable interpolation between time grid points. Then we prove strong convergence with order 
7 = 1/2 under some regularity assumptions on the coefficients of the SDE. In section 3, we 
propose a modified Ninomiya-Victoir scheme, which may be strongly coupled with order 1 to 
the Giles-Szpruch scheme. This result allows us to derive an antithetic version of the Ninomiya- 
Victoir scheme and combine the ideas of Giles-Szpruch and Debrabant-Rossler by building the 
multilevel Monte Carlo estimator with the Giles-Szpruch scheme from level 0 to level L — 1 
and the coupling between the Ninomiya-Victoir scheme and the Giles-Szpruch scheme at the 
last level L. The efficiency of this estimator is confirmed in section 4, where we present and 
comment, in details, numerical experiments carried out on the Clark-Cameron SDE and Heston 
SDE as in [ 6 ]. 

2 Strong convergence of the Ninomiya-Victoir scheme 


We begin this section by introducing some notations which will be used throughout this paper. 
To discretize (1.1) we consider a uniform grid with time step h = T/N where N gW and we 
denote: 


• {tk)kelo-N} ~ subdivision of [0,T] with equal time step h, 

• fs the last time discretization before s G [0, T], ie = tfc if s G (tfc, tfc+i], and for s = to = 0, 
we set To = to = 0, 

• fs the hrst time discretization after s G [0,T], ie fg = tfc+i if s G {tk,tk+i], and for 
s = to = 0, we set tq = to = 0, 

• Vj G {1,..., d} , Vs G [0, T], such that tk < s < tk+i, AWi = Wi — 

• Vs G [0, T], such that t^ < s < ffc+ij As = s — t^, 

• T] = (iji,... a sequence of independent, identically distributed Rademacher random 
variables independent of W, 
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• By a slight abuse of notation, we set = rjk+i if s G {tk,tk+i], 

• Vx G M+, \x\ denotes the unique n G N* satisfying n — 1 < x < n, 

• Vx G M+, [xj denotes the unique n G N* satisfying n < x < n + 1. 


Let V : M” —)• Lipschitz continuous and consider the ordinary differential equation in M"': 

/ ^ = (2 1) 
lx(0)=xo. 

The solution of (2.1) at time t, t G M is denoted by 

x{t) = exp{tV)xo, 

and the integral form of (2.1) is given by 

x(t) = exp(tB)xo = xq + / V{x{s))ds = x+ / B (exp(sl/)xo) ds. 

Jo Jo 

We recall that in (1.1), each coordinate i G {1,..., n} evolves according to the following stochas¬ 
tic differential equation 

d 

dXi = b\Xt)dt + (T^^{Xt)dWi. 
i=i 

Then, assuming regularity for the diffusion coefficients, one can write (1.1) in Stratonovich 
form 


dXt = a^{Xt)dt + Y. o dwi 

i=i 

Xq = x 

d 

where = b — ^ Y da^a^ and da^ is the Jacobian matrix of defined as follows 
i=i 

da^ = 

Now, we present the Ninomiya-Victoir scheme introduced in [8]. 

• Starting point: X^’^ = x. 

• For A: G {0 ..., — 1}, if r/fc+i = 1: 

.NV,'r) _ f h. 


( 2 . 2 ) 


Kl. = • • • «"P ) «-p ( 

and if -qk+i = -1: 


tk ’ 


= «"P { ) ^"P ) • • • ^-P () exp I ) W, 




(2.3) 


(2.4) 


The Stratonovich form is preferred when we use the Ninomiya-Victoir scheme since the Stratonovich 
drift appears in the dehnition of the scheme. Moreover, using Ito’s formula, one has 

Vt, s G K+, s < t, exp [{Wl - Wj) v) y = y + j V (exp {{Wj, - Wj) V) y) o dW^ (2.5) 
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Then, rewriting (2.3) and (2.4), one obtains 


^ + l r^k-\-l 1 / / \ \ 

<+? = + E / ° + I 2 (^° (2-6) 

where, for s G (4,4+i], 

X"-” = exp , (2.7) 

for s G {tk,tk+i], j G {1,.. .,d}, 

Xi” = exp (XWia>) , (2.8) 

for s G (4,4+i], 

= exp ■ (2.9) 

Denoting = ^tu+i^ ~ ^ets an expression similar to (2.8) for j G {0,(i+ 1} 

and s G (4,4+i] 

xi" = exp (^.x») ■ (2.10) 

Then, one can observe that the Ninomiya-Victoir scheme is obtained by replacing the exact 
solution X by one of the intermediate processes in the Stratonovich formulation (2.2) of 
the SDE (1.1). 

Remark 2.1 The stochastic processes ^ j G {!,... d + l}, are not adapted to the 

natural filtration Ft = a{Ws-,s < t) of the Brownian motion. To get around this problem, we 
work with the following filtration F^ = a (^Wi, s <t) \/ a [W^, s < T) , Vj G {1,..., d}. Then, 

for j G {1,..., d}, by independence, is a F^ Brownian motion, and X^’"^ is adapted to the 
filtration FT This ensures that each stochastic integral is well defined. 


In order to study the strong convergence, we have to build an interpolated scheme. Let 
be the following Ito process 


XNV,,) 


/te[o,r] 


dX 


NV,-n _ 


0 f yd+l,ri 


.t = E ^Hxn o dWi + 1 (a^ X)^’^ +a^X, 
^=1 ^ ^ ' 

Xo^^'^ = X. 


dt 


( 2 . 11 ) 


Using (2.6) and forward induction, one can show that is an interpolation of the 


0<t<T 


Ninomiya-Victoir scheme ^ The Ito decomposition of ^ ^ is given by 


te[o,T] 


dXr^^ = E ^KxndWf + i E da^a^ {Xfi^]dt+^^ [a^ (X^) + (X 


XE" = X. 


1=1 


1=1 


d+l,r) 


dt 


( 2 . 12 ) 
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Remark 2.2 A natural and adapted interpolation for this scheme could be 


where 

/i_i (to, .. .,td+i;x) = exp (ioo-°) exp (ticr^) .. .exp (^tdcr‘^^ exp (trf+ia°) x (2.14) 

and 

hi (to,td+i]x) = exp (iocr°) exp ... exp (iicr^) exp {td+icr^) x. (2.15) 

In both cases AWt = (AW/,..., AVk/). In order to obtain the ltd decomposition of , 

we have to apply the ltd formula. To do so, we have to compute the derivatives of /i^. In the 
general case, the computation of derivatives of this function is quite complicated. That is why 
we will not focus on this interpolation. 


(2.13) 




— ,AW,— 

2 ! 1 ) 2 ’ Tt 


2.1 Strong convergence 


We recall that {W^, M"), Vj G {1,..., d}, and we assume that the vector fields, cr-^ , Vj G 

{0,...,d}, and da^a^,\/j G {l,...,d}, are Lipschitz continuous functions. Obviously, b is also 

d 

Lipschitz continuous, since b = ^ da^aT Let L G denote their common Lipschitz 

i=i 


constant: 


Vj G{0,...,d},V3:,?/GM”,||o-^(a:)-o-^(y)|| <L\\x-y\\, 

Mj G {1,..., d} , Vx, y G M"", a\x) — da^a^{y)^ < L\\x — y\\ 


where the Euclidean norm is denoted by 


Theorem 2.3 Let p G [l,+oo). Under the previous Lipschitz assumption, there exists a deter¬ 
ministic constant Cmv £ 1^1 such that 


VAgN*, e 

sup 

Xt - xf^’^ 

2p 

V 


t<T 





< Cnv ^1 + h^. 


Of course, this result implies that 


VAg N*, 


E 


sup 

t<T 


Xt-X, 


NV,ri 


2p 


< Cj\fv (^1 + 


Obviously, ^ and h depend on N, but in order to keep the notations simple, the 

dependence on N is not made explicit. The following proposition will be used to prove the 
theorem. 


Proposition 2.4 Letp > 1, T = (Lt)o<t<h ^6 the solution of the following n-dimensional SDE, 
driven by a d-dimensional Brownian motion, until t = h 


dYs = a{Ys)ds + P{Ys)dWs 
Yq independent of (VLi)^gjQ such that E 



< + 00 . 


Assume that a and /3 are Lipschitz continuous functions, then: 3Co G Vt, s G [0, h], s < t, 
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(^) 

(ii) 


E 




< E 


E 


\Zt — Zs\ 


\2p 


< Cq (1 + e 


l + ||Zo||^^ exp(C'o/i). 

) (i - «)’ 


\Znfp 


If f3 = 0, we have a better result: 


E 


\Zt-ZsfP 


<C'o(l+E[||ZofP])(t-s)^P. 


(2.16) 

(2.17) 

(2.18) 


The constant Cq only depends on ||a(0)|| 
tions a and /3. 


, T, p, and the Lipschitz constants of the func- 


All these results are well known (see [11] for example). 


2.2 Intermediate results 

By using the previous proposition, one can show that the scheme has uniformly bounded mo¬ 
ments. 

Lemma 2.5 Vp > 1,3(71 £ Wf^t G [0,r],V77 G N*,Vj G {0,..., d-I- 1}, 





2p 


E 

1 + 



V 


< exp((7ift) (^1 -t . 


Proof : Let p > 1 and t G [0,T], then 3A; G {0,... , iV — 1} such that t^ <t < tfc+i- For j = 0, 






is the solution of the following ODE 


dZg = ^a^{Zs)ds 


_ yNV,ri 


Zo. = X 


The independence between p and W combined with (2.16) ensures that: 


E 


1 + 


2p 

V 

< E 

1 + 

_ 









2p 


exp ( -Coh 




1 + 


X 


NV,r] 

^k 


2p 




1 + 


X}’^ 

tfe + 1 


2p 


exp ( -Coh 


(2.19) 


Similarly, for 1 < j < d: ( 


is the solution of the following SDE: 
dZs = Ida^a^ (Z,) ds + a^{Zs)dWi 


Ztk = X 


tfc + 1 
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Using the same argument, one gets: 


E 


1 + 




2p 




1 + 


Xj-i,v 
tfe + 1 


2p 


+ l{%+i=-i}IE 


1 + 


Xj+i,v 

tk+l 


2p 


exp (Cq/i) 


( 2 . 20 ) 


Obviously, for j = d + 1, one has a similar result: 


E 


1 + 


xr 


2p 




1 + 


xd,^ 

C/c + l 


2p 




1 + 


X 


NV,p 


2p 


exp ( -Coh 


( 2 . 21 ) 


The global Lipschitz constant L is the same for all vector fields, therefore, the same constant 
Cq is involved in the three inequalities. In both ODEs, the vector field is multiplied by 
1/2, it is equivalent to integrate the equation until h/2 and simply remove the multiplicative 
factor 1/2. That is why one gets a factor 1/2 in both inequalities (2.19) and (2.21). Since for 
all A: e {0,..., N}, X^^’^ = one can use forward induction on 

k combined with forward induction (respectively backward) on j G {0,..., d + 1} if r]k+i = 1 
(respectively = — 1) to get: 


E 


1 + 


xr 


2p 


< exp (Citk+i) 1 + ||x 


|2p 


where Ci = (d + 1) Cq- 


The following lemma is a direct application of Proposition 2.4, together with Lemma 2.5. 
Lemma 2.6 Vp > 1,3(72 G Vf G [0, T], VX G N*, Vj G {1,... ,d}, 


E 


- Xf^ ^{vt=i} - Xf, ^{pt=-i} 


2p 


< (72 (l + \\xfp^ hP, 


and for j G {0, d + 1}, 


E 




2p 


<02(1 + lixfp) rp, 


where by convention X-^’^ = xf^'^ = X^^’^. 


Proof : Let p > 1, t G [0,T] and j G {1,... ,d}. Thanks to (2.17) in Proposition 2.4 we have: 


E 




2p 


< (7o 1 + E 




2p 


Since 


1+E 




2p 


l{^t=i}E 


1 + 




2p 


+l{r?t=-l}IE 


1 + 


X: 


j+^,v 


2p 
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combining this estimation with lemma 2.5 we get 


E 




2p 


< C'oexp(C'ift) (^1 + 

< Coexp(CiT) (^1 + hP 


Applying a similar argument, using (2.18) from Proposition 2.4, we get the same result for 
and We conclude by setting C 2 = Cq exp (CiT). I 


The following lemma deals with the estimation of the difference between the scheme and 

the intermediate process for j G {0,..., d + 1}. 


Lemma 2.7 Vp > 1,3(73 £ G [0,T],VA G N*,Vj G {0,..., d + 1}, 


E 


A/’'' - A 

t n 


NV,r] 


2p 


<(73(l + ||a:f^)dP. 


Proof : Let p > 1, t G [0, T] and j G {1,..., d + 1}. Using telescopic summation and convexity 
inequality, we get 


i n 


2p 


< (d +! 




2p 


+ E 

ptm<ptj 


^ft - ^ft Mvt=i} - Xft Hvt=-i} 


2p 


Taking the conditional expectation, and using Lemma 2.6, we obtain 

2p 


E 


t ^Tt 
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< 


(d + (d + 2TP) C 2 (1 + hP = (73 (1 + hP 


with (73 = (d + 2fP-^ (d + 2TP) C 2 . 


2.3 Proof of the strong convergence 


Proof : Let p G [l,+oo), t G [0,T] and s G [0,t]. Subtracting (2.12) from (1.1), we can 
evaluate the difference between the exact solution and the scheme 
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Using a convexity inequality and taking the conditional expectation of the supremum, we get: 


E 


sup ||X, - 


2p 


s<t 


d 


1 


d+l 


< (2 (d + 1))^^ ^ ^ 

.i=i j=o 


( 2 . 22 ) 


where 


/o = E 


Id+l — IE 

and for j S {1,..., d} 


Ej — E 


sup 

S<t 


sup 

s<t 


na\Xu)-a^ {X°u^))du 
Jo 


sup 

s<t 


I,=E 


sup 

s<t 


r {a^{X^) - a^iXi’^)) dWi 
Jo 

[ {da^a^ (Xu) - da^a^ (^d’’')) du 

Jo 


V 


Let us focus on Ej and dj, for j E {1,... ,d}. The independence between W and r] permits to 
apply the Burkholder-Davis-Gundy inequality to obtain 


Ej<KE 


/ ft o \ ^ 


ft 



/ hHx..) - aUxim du) 

V 

< KTP-^ / E 

\\a^{Xu)-a^(Xn\\"" 

V 

\Jo J 


Jo 




du 


where K is the constant that appears in the Burkholder-Davis-Gundy inequality. By the Lips- 
chitz assumption 


Ej < KTP- 


Applying a convexity inequality, we obtain 


ij < r^p- 


Again, by the Lipschitz assumption, we also get 

rt 


[ E 

\\x - 

II 

_1 

Jo 




du. 


(2.23) 


f E 

||5ctV^' (X,) - (Xf’^) 

_1 

Jo 




ds. 


Ij < t^p-^L'^p [ 
Jo 


E 


X - X^’P\ 


|2p 


du. 


(2.24) 


Using the same approach, we get a similar result for Iq and Id+i- Gombining (2.23) and (2.23), 
together with (2.22), we obtain 


E 


d+l 


E 


X - X^'P\ 


|2p 


du 


sup\\Xu-X^^’P\\‘^^ T] <a'^ [ 
s<t j=0 

where a = (2 (d + 1 ))‘^p~^ L‘^p (^KTP~^ + . Now, we look at 

{0,..., d-|-1} and u £ [0, t]. Introducing the solution X and the Ninomiya-Victoir scheme 
X^V^Tj using a convexity inequality we get 


Xu - X^I^ 


(2.25) 
Let j e 


E 


X — 


|2p 


< 3^P~^ E 


\Xu-Xf„fP + 


NV,ri 


Xi- -X. 

Tu Tu 


2p 


+ 


xNV,v _ xlv 


2p 


10 









































































Then, using the estimation (2.17) from Proposition 2.4 


E 


\Y — Y- W^P 


< 


Co (^1 + (u - TuY < Co (^1 + ||a:||^^^ hP 


and from Lemma 2.7 


E 


- X, 


NV,v 


2p 


< Cs (l + llxf?’) hP. 


Moreover 

We finally get 
E 

where 

and 


E 


-X 

Tti Tu 


NV,T] 


2p 


< E 


sup ||X^ - xY^’Y 


2p 


V<U 



V 

< /3 / E 


V 

du + j(l+ X 

s<t 


do 

v<u 


\ / 


13 = {d + 2) 


a 


7 = /3T(Co + C3). 

Before applying Gronwall’s lemma, let us remark, by (2.25),(2.16) and Lemma 2.5, that E 


sup 

y ^NV,r) 

■^s 

_ s<t 



2p 


is finite. 

We conclude thanks to Gronwall’s lemma 


E 

sup 

X, - 

2p 

V 


t<T 





< exp (/3T) 7^1 + hP. 


We conclude this section with a lemma which will be useful for the next section. 


Lemma 2.8 Let F G (M”,M”') and assume that its first and second order derivatives have 
a polynomial growth. Under the assumptions of the Theorem 2.3 we have the following result: 
VpG [1,+oo), 3C4 GM;,Vi G {0,..., d + 1} , VX G N*, 




pt / Ti TT r \ 

2p 

■ 

E 

sup 

/ F {Xi-") - F (x3X) 


V 


t<T 

Jo ^ ^ 




< C^h^P. 


Proof : Let j G {0,..., d + 1} , z G {1,..., n}, and t G [0, T]. Using the integration by parts 
formula 


[ (f- {X}”) - F‘ (x"*'’”)) du = r (t A f. - s) d (F‘ {X}”)) + r Yi ‘‘ T (L"*’)) ■ 


Psm<psj 
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Then, using the chain rule for m G {0, d + 1}, we get 

d (F* . VF* ds. 

Applying Ito’s formula for m G {1,..., d}, we obtain 

d (F* (X™’^) . VF* (X™’^) + ifr (a”* (a”*)* (X™’^) V^F* (X™’^))^ ds 

+ . VF* (Xf^’^) dW^. 

In both cases, combining a convexity inequality, the Burkholder-Davis-Gundy inequality, the 
Holder inequality, the Lipschitz assumption on cr"*, da'^a™' for m G {0,... , d}, the polynomial 
growth assnmption for the first and second order derivatives of F, and t/\fs — s < h, Vs G [0, ft], 
we get two constants 7 G Ml and g G N*, independent of N, such that 


E 


sup 

t<T 


F* (x^’^) - F* 


ds 


2p 


rf+l „T 


< 


f 

m=oJO 


E 


1 + \\X^’^ 


\2q 


ds. 


We conclude by using Lemma 2.5 and taking the Euclidean norm. 


3 Coupling with Giles-Szpruch scheme 


In [ 6 ], Giles and Szpruch proposed a modified Milstein scheme defined as follows 


i (^£h (‘*^+1 - *0+E(-^<1) 


i=i 




j,m=l 


(3.1) 


vGS _ 
^to - X. 


rtk+1 

In comparison with the Milstein scheme, the terms involving the Levy areas / AW^dW^ — 


C^k + l 


’ ik 


^k 


AHVdlT^ have been removed. According to Lemma 4.2 in [ 6 ], the strong order of con¬ 


vergence is 7 = 1 / 2 . 


Lemma 3.1 Assume that b,a^ G C^(M”,M”),Vj G {!,...,d}, with bounded first and second 
order derivatives, and that da^a'^ ,Vj, m G {1,..., d}, has bounded first order derivatives. Then: 


3Cgs g m;,vx g m;. 


E 


max 

k&{0,...,N} 




vGSw'^P 

^tk II 


< CcshP- 


(3.2) 


Giles and Szpruch also proposed an antithetic version of their scheme based on the swapping of 
each successive pair of the Brownian increments in the scheme. With regards to the multilevel 
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Monte Carlo estimator, Giles and Szpruch use the arithmetic average of the scheme (3.1) and 
its antithetic version on the fine grids, at each level I G {1,... ,L} as follows 



The swapping of each successive pair of Brownian increments provides a strong convergence of 
order 1 between the schemes used on the coarse and fine grids, and so Giles and Szpruch obtained 
the convergence rate /3 = 2 of the variance V — / (^X^ 

when the payoff / is smooth. In this way, using this multilevel Monte Carlo estimator leads to the 
computational complexity O (e~^) for the mean-square-root error e. To use the Ninomiya-Victoir 
scheme either at each level or only at the finest level of a multilevel Monte Carlo estimator, we 
study in this section the coupling between the Ninomiya-Victoir and Giles-Szpruch schemes. 
To keep /3 = 2, we suggest to compare the Giles-Szpruch scheme with the following modified 
Ninomiya-Victoir scheme 

xNV,v ^ 1 . (3.3) 

To be consistent with the interpolation of the Ninomiya-Victoir scheme, we define the interpo¬ 
lation of the scheme between the grid points as follows 





Theorem 3.2 We assume that b £ (M”;M”) with bounded first and second order derivatives, 

and G (M”;]R"') ,Vj G {1,... ,d} , with bounded first and second order derivatives and with 
polynomially growing third order derivatives, and that da^a^, \/j,m G d}, has bounded 

first order derivatives. Then: 


3C G M;,W G N*, E 

sup 

XNV,v _ x^S 

2p 

V 


t<T 





< Ch^P. 


Proof : We denote by L the common Lipschitz constant of b, and m G {1,... , d}. 

We also denote by M the global bound on first and second derivatives of b and afiyj G {1, ..., d}. 
Let t G [0,T] and s G [0,t]. Writing X^^'P in integral form, we get 



(Xi-’>))iWi + Y. -.{Sc 


{XfiP) + da^a^ {X^y-^)) du 
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Then using \h — da^a^ — = 0, we get 

i=i 


]+b{X 

Tu ) \ Tu 


xr -=-+E K"'-’)) -nxi + /K" ( 

j _ 1 0 </ 0 

+ E /’ 5 ("' (Xt”) - ai + ai {Xj->) - <,= {x^X")) dWl 

j=l JO 

+ E I i^u-^) - (<'"’■")) 

i=i 

+ i (<r» {XX - ‘t” {x?X) + (■?:’-”) - du 

+ /■ j (<^“ (i?'-") - <^“ (xX”) + ( 

0 

Subtracting (3.4) and using a convexity inequality, we obtain 
E sup||xf^’^ 


NV-tj 


du 


du 


' Xf ) ) du. 


d d d-\-l 

< S^f-' (d + I)"*-* I ^ £j + ^ R, 

j=l j=0 j=0 


(3.5) 


where 



2p 


2p 


Jid+l = E 


sup 

S<t 


i (<r» (xX-”) - <r° (xfX) + (XX-”) - <r° «''’-’)) du 


2p 


and for j G {1,..., d}, 


Ij = E 


sup 

S<t 


i {Xi”) - a’ [xX”) + ui {XX) - <r’ {xX'")) d^i 



da=a= (XlA AM^- + - 53 a.rV’" (Xg®) AW^ dWl 


m=l 

m^j 


2p 


Ej — E 


sup 

S<t 


Rj = E 


sup 

S<t 


1 [ad [xX") + (Vg''-’)) - (Xff)) dWi 

i (^da^aJ (Xi’^) - da^aJ + da^du 


2p 


Step 1: estimation of Ej, for j G {0,..., d}. 

Let us start with the estimation of Ej, for j G {0, ... ,d}. We set for = b and F^ = 
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for j G {1,... ,d}. Combining the Burkholder-Davis-Gundy inequality and a convexity inequal¬ 
ity, we get 

2p 


E, < ma^{T^P-\KTP-^} E 1 (f^' {xgS) 


V 


du 


where K is the constant that appears in the Burkholder-Davis-Gundy inequality. For i G 
n}, denoting Yu = ^ [f^^ performing a 

second order Taylor series expansion, we obtain 

u = {K'''") -^ {YY - YF")' (ft) + (fl)) {xY” 

where and are points between and Then, we easily get 

||2d ^ Wr, 2p ..An /»1 ._An/_n '^P' 

|F„|| ^ <ai 


- X 


NV-ri 


( 

j^NV.p _ ^GS 

2p 

-b 

XNV,r) _ x^Y-P 

V 

Yu '^u 


Yu Yu 


where ai = 2^^ ^ -|- • 

Fj < aimax{r2p-^FrP-i} ( [ 


Thus 


E 


sup 

V<U 


\XNV,v _ j^GSfP 


du+ E 

Jo 


X^^’^ - X 

Yu Yu 


NV-t) 


A.p 


du 


Introducing the solution X at time f„ and using a convexity inequality, we obtain 


E 


X^v,v _ x^Y-v 


4p 


< 2^P 


-1 


E 




4p 


-hE 




NV-ri 


4p 


Thanks to Theorem 2.3, we deduce that 


E 


XNV,r] _ x^Y-V 


Ap 


< 2‘^pCnv (l + h^P. 


It follows that 


F, < /3i 


'i: 


E 


sup 

V<U 


I vW,p vGSw'^P 
A,, — A.,, 


du -|- /i^P 


(3.6) 


where/3i = ai max {r^P ^,KTP max |i, 2^PC'Ary -|-||x||^^^ |. 

Step 2: estimation of Rj, for j G {0,..., d}. 

Turning to the estimation of Rj, for j G {0,..., d}, from Lemma 2.8 we get a constant /32 S 
such that 

Rj < hh^P. (3.7) 

Step 3: estimation of Ij, for j G {1,..., d}. 


It remains to estimate Ij, for j G {1,... ,d}. Using the Burkholder-Davis-Gundy and convexity 
inequalities, we get 


I, < 


E 


a' {Xi!") - a’ [X^X) + Yt") - (v"'’-’) 


2daGi (Xg^) XWi - Y, (Vg®) AW, 


m=l 

m^j 


2p 


ds. 
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Introducing 


U/J = I xhJ-V 

^ u ^ u ' ^ u 


where 


- da^a^ da^a^ (3-8) 


ri„m<riuj 


and 


rium<riuj 


= da^a^ AWi + da^a^ XW^ - 2da=a^ {xg^) XW^ + [xgg^'^) AVT, 

d 

+ ^ 9cjV™ (xggg'-A XWg - ^ {xfS) XWi 


m 

fu 


/■m 

fu 


r)um>r}u3 


we obtain 


h < -KTP 

J — rs 


-1 


m=l 

mj^j 


E 


(3.9) 


l^i 11^?’ + ||$i \gP 


du 


Step 3.1: estimation of E 




2p 


, for j e {1,.. .,d}. 


Applying Ito’s formula in (3.8) to compute [xt^^ — cr^ , we get 


= 


da^a^ {Xi’'^) - da^a^ (xg'^’"’)) dW^ + {da^a'^ (A™’’?) - da^a^ [xgg^'^'Yj dW, 


m 

V 


VurrKriuj 


+ l r {Xi’^) dv ^ H OF^’^a^ (A)^’'’) dv 


where F^’'^ = da^a'^ for j, m G {1,..., d}. Note that the term ^ X) / 9F-^’”*(T™ (X^’’^') dv 

rjurrKrjuj Jx 

is equal to the sum of the drift contribution and the ltd correction due to the dynamics of . 
The assumptions on and Lemma 2.5, ensure that 

• Vj, m € {1,..., d} , is Lipschitz continuous. 

• Vj, m G {1,..., d} , (A™’^) has uniformly bounded moments. 


Using Lemma 2.7, the Burkholder-Davis-Gundy and a convexity inequalities, we obtain a con¬ 
stant 73 G M;)_ such that 

|2p 


E 


U ' 


r] 




Obviously, we have the same inequality for 

2p 


Step 3.2: estimation of E 




, for j G {1,... ,d}. 
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By the Lipschitz assumption 


< {d+lfP~^L^P\ 


j^NV,r] _ j^GS 


2p 


\AWif’’ + 


NV-n _ j^GS 


+ E 

ijuirKriuj 

By independence 


^NV,t] _ j^GS 

Tu '^u 


2p 




+ E 

Pum>puj 


X. 


j^N'V^—r] _ 


2p 


\^wi\ 


2p 


2P , 




(3.10) 


E 


j^NV,v _ j^GS 




= E 


j^NV,v _ j^GS 


2p 


E 


Aiyi 


< 2^P"^E 


AW^‘ 


E 


- Xf 

Tu 


2p 


+ E 


Iv yGSw'^p 

\^Tu II 


Then, using Theorem 2.3 and Lemma 3.1, we get 


E 


j^NV,v _ j^GS 

Tu '^u 




V 


< 2^P E (Cnv ( 


1 + ||x||^^ ) + 


C'gs) 


where G is a normal random variable. Using the same approach, we get the same result for the 
other terms in the in the right-hand side of (3.10). Thus, we deduce that there exists a constant 
as e Ml such that: 


E 




) 112P 


V 


< ash'^P. 


(3.11) 


Combining our different inequalities, we obtain 

Ij < 

where f3s = ^KTP (as + 2 ^P 73 ) . 

Step 4: conclusion 


Finally, by combining (3.6), (3.7), (3.11), together with (3.5), we complete the proof using 
Gronwall’s lemma 


E 


t<T 


V 


< Ch^P 


where C = 3^^ ^ {d + VfP ^ {dfii + {d + l)/32 + (d + 2) f5s) exp ^3^^ ^ (d + 1)^^ ^ d/3iT^ . 


4 Multilevel methods for SDEs 


In this section, we are interested in the computation, by Monte Carlo methods, of the expectation 
y = E [/ (At)]) where X = (At)jgjQ is the solution of the stochastic differential equation (1.1) 
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and / : 


f-)- 


a given function such that E 


fiXT? 


is finite. We will focus on minimizing 


the computational complexity subject to a given target error e. To measure the accuracy of an 
estimator Y, we will consider the root mean square error 


/ ^ \ 1 


2" 

[Y,Y] = E2 

Y-Y 



4.1 Multilevel Monte Carlo 


The multilevel Monte Carlo method, introduced by Giles in [5], consists in combining multiple 
levels of discretization, using a geometric sequence of time steps hi = T/2^ for example. Denoting 
by a numerical scheme, with time step T/N, the main idea of this technique is to use the 
following telescopic summation to control the bias 

= E [/ (X^)] + J^E [/ ■ 

1=1 

Then, a generalized multilevel Monte Carlo estimator is built as follows 

L Ml 

^MLMC = (4.1) 

/=0 ^ ■'* k=l 

where (^fc)o<Ki, i<k<Mi independent random variables such that for, a given discretization 
level Z € {0,..., L}, the sequence {^1) is identically distributed and satisfies 

E [ZO] = E [/ (X^)] (4.2) 


E 




and 


V/€ {1,...,L},E 




E 


/ U: 



(4.3) 


Assume that, for a given discretization level / G {0,..., L}, the computational cost of simulating 
one sample is C'A;2*, where C G M+ is a constant, depending only on the discretization scheme 
and A; G Q+ is a weight, depending only on 1. The computational complexity of Ymlmc, denoted 
by Cmlmc, is given by 

L 

Cmlmc = C E MiXa^. (4.4) 

1=0 

The natural choice for £ {0, ... ,L} considered in [5] is 


Z^ = f (X^) 


(4.5) 


Z^ = f (Xl') - / , VZ G {1,..., L} . (4.6) 

For this canonical choice, it is natural to take Aq = 1 and A; = 3/2, V/ G {1,..., L}. According to 
Theorem 3.1 in [5] the optimal complexity depends on the order a of weak convergence 

of the scheme and the order (3 of convergence to 0 of the variance of ZK Here, we recall this 
complexity theorem. 
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Theorem 4.1 Assume that 


E 


f Xj 




and 


^ ( ^0 2 / 5 « 


(4.7) 

(4.8) 


for some constants ci G M* and C 2 G Ml independent ofl. Then, by choosing: 


L* = 


log2(^ 


a 


(4.9) 


and 


wg{o,...,l*},m; = 




j=0 


(4.10) 


we get an optimal computational complexity: 

^MLMC = O {e if /3 > I 


MLMC — Ole ^ ( log ( - 


MLMC 


= 0 e 


- 2 + 


S-l 



if I3 = l 


(4.11) 


if 13 < I 


with RMSE ( Ymlmc, Y ) bounded by e. 


To obtain the estimation (4.8), the key point is that the simulation of X'^‘ and ^ comes 
from the same Brownian path. We easily bound the variance convergence rate from below using 
the strong convergence rate 7 of the numerical scheme, since in general, /3 > 27 for a smooth 
payoff. To attain 7 = 1 , one has, in general, to simulate iterated Brownian integrals involving 
Levy areas, for which there is no known efficient method. To get around this difficulty, Giles and 
Szpruch introduced a Milstein scheme without Levy areas and its antithetic version by swapping 
the Brownian increments. In a multilevel Monte Carlo method, using the arithmetic average 
of the modified Milstein scheme and its antithetic version in the hnest grid, and the modified 
Milstein scheme in the coarsest grid leads to /3 = 2. By this way, Giles and Szpruch managed to 
improve the variance convergence rate without simulating the Levy areas. To be precise, they 
choose as follows 


Z^GS= /(^' 


.GS,1 


7* — 

^GS — 




x: 


GSY 


+ / ( 


-f{X^^’^‘~"),'^lG{l,...,L}. 


(4.12) 

(4.13) 


Here, is the Giles and Szpruch scheme defined by (3.1) using a grid with time step 

hi = r/2' and is an antithetic discretization defined by swapping each successive pair 

of Brownian increments in the scheme. To be more precise, we dehne two grids, a coarse grid 
with time step /i/_i and a fine grid with time step hi. The discretization times (ifc)o<fc< 2 *-i 

^ j j are defined by tfc = A:/i/_i,V/c G { 0 ,..., 2^“^} , and = (/c + ^) G 
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r 7—1 '1 ( 

|0,..., 2 — Then, on the coarsest grid, ) k { ' defined inductively by 


= X and 


^GS,2'-i , , / ^G5,2'-i 
+01^. 






j,m=l 


.GS,2‘ 


i=i 




where - Wt^. Similarly, on the finest grid, 


(4.14) 


is defined 


GS 2 } 

inductively by ’ = x and 


K'S = + b (^tf’"') hi+E ’"') AIT/’^ 




G5,2* 

tk+i 


i=i 


+4 X: A, 


j,m=l 


.GS,2* 

'-4 


AM/" AH',™/ - 1, 


Tc+J 




i=i 


k +12 

r)l 

)' 

"fe+^ 


+hE^»</"” (am/" 


(4.15) 


where AWf , = Wt , — Wu, AW/ , = w/ , — w/ , . The antithetic scheme is defined 

T+l T+l ''=’ ifc+l ifc+l T+1 

by the same iterative equations, except that the Brownian increment AWf and AW/ are 

swapped 


X 


GS,2‘ 

^k+^ 




^GS,2' 
‘^4 


J=1 


^GS,2* 

“^4 


AM/". 


+i E (If •^') (AM/f.AM/f-; - 1, 

j,m=l \ / \ / 


x: 


GS,2' _ yGS,2^ 

tk + l 


rGS,2‘ 


k+i 


=Kit +' ( 


j=l 


k+4 


k+4 


(4.16) 


Theorem 4.10, Lemma 2.2 and Lemma 4.6 in [6] ensure that (3 = 2 under some regularity 
assumptions on / and the coefficients of the SDE. 


Theorem 4.2 Assume that f £ C^(M”,]R) with bounded first and second order derivatives, 
b,a^ £ C'^ {EA,EA) ,\/j £ {l,...,(i}, with bounded first and second order derivatives, and that 
da^a'^, yj,m £ {1,... ,(i}, has bounded first order derivatives. Then: 


Vp > 1,3c £ m;,V/ £ N*, E 



2p 



where Zq^ is defined by (4.13). 
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To account for the use of three schemes in the levels I G {1,..., L*} instead of one in level 0 
we choose Aq = 1 and \i = 5/2, VZ G {1,... ,L*} . Then, the multilevel Monte Carlo estimator 

^MLMC ~ '^^GSi where L* and M* are given by (4.9) and (4.10), respectively, achieves 
1=0 ‘ 

a complexity O (e~^)- In [2], Debrabant Rossler improved the multilevel Monte Carlo method 
by using, in the last level L, a scheme with high order of weak convergence. Although this 
modified method attains the same complexity, it reduces the computation time by reducing the 
bias. We can follow this idea using the Ninomiya-Victoir scheme at the last level L, thereby 
taking advantage of its order 2 of weak convergence. More precisely, we propose to choose 


Z^GS = f[X^ 


Zhs = l^ 


GS,2‘ 


^GS-NV 


l(/(- 


X 


NV,2^,rj 


.GS,2‘ 


+ / X 




X‘) 

(4.17) 


(4.18) 




(4.19) 


Here, (respectively ~^) is the antithetic discretization of the Ninomiya-Victoir 

scheme (respectively ’~^)) obtained by swapping each successive pair of Brow¬ 

nian increments. Theorem 3.2 ensures that (4.8) the order of convergence of the variance at the 
last level L is 2. 


Proposition 4.3 We assume that f £ C'^ (M”',M) with bounded first and second order deriva¬ 
tives, b G C^(M"’,M”) with bounded first and second order derivatives, G C^(M”,M"’), Vj G 
{1,... ,d}, with bounded first and second order derivatives and with polynomially growing third 
order derivatives, and that da^a^, Vj, m G {1,..., d}, has bounded first order derivatives. Then: 


Vp> 1,3c G M;,VZ G N*, E 


^GS-NV 


2p 



where ZQg_j^y is defined by (4.19). 


Proof : Let p > 1, introducing ^ ( / ( 
we get 




GS-NV 


2p 

- 22P 


5 (/(- 


NV,2\p 




and using a convexity inequality. 




X 


NV,2\p 


+ / X 


:nv,2\-p 




2p' 




-1 


yl 

^GS 


2p 


TT / vNV,2h'n v^V,2‘,-r] ^GS,2* 

However, [Xy jX.p 


and (X 


NV,2hp yNV,2h-v yGS,2'' 


I Y 

, A, 


, Xy ’ I have exactly the same 


distribution. Then, by taking the expectation we obtain 


E 


■'GS-NV 


2p 


o2p-l 

< V-plE 

- 22p-l 


X 


NV,2\rt 


+ / X 


rNVX-p 


- f (X^^’2' 


2p 


-T 3^P"^E 


^GS 


2p 


Denoting Xj 


YNV,2‘,ri ^ 1 fj^NV,2‘,p ^ yNV,2\-ri 


XI 


and performing a second order Taylor expansion 


as in Lemma 2.2 in [6], we get a constant C G which only depends on / and p, such that 


E 


X 


GS-NV 


2p 


< C E 


^NV,2hp vGS,2'' 

yXrp XYrp 


2p 


-fe 


vNV,2\p vNV,2h-p 

u\. rp ./J. rp 


4p 


-h E 


X 


GS 


2p 
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Introducing the exact solution X at time T in E 


E 


X 


NV,2\r] yNV,2‘,- 


4p 


< ( E 


E 

vNV,2\p 

p 

-X 

:Nv,2‘,p _ 

4p 


4p 


, we get 


+ E 


Xt - X^^’2 


4p 


Since ’^,Xr^ and ’ ^,Xt^ have the same distribution, we deduce that 


E 


yNV,2\r) vNV,2‘,-ri 

y\. rp ./J. rp 


4p 


< 2^PE 


X 


NV,2‘,ri 


-Xt 


4p 


Hence: 


E 




GS-NV 


2p 


< 2^PC ( E 


yNV,2'',ti yGS,2^ 

v\. pi p 


2p 


+ E 




Ap 



2p 


-7E 

Zhs 



Then we conclude using Theorems 2.3, 3.2 and 4.2. 


Exploiting the telescoping summation, one can change the constraint (4.3) on the last level L 
and assume: 

(4.20) 


r2^ 

'-T 


E [Z^] =E\f(x^)-f (X 

Here X is an other scheme, and to be consistent, (4.7) becomes 


E 


/ U: 




L-l ^1,1. T l. 

Then we propose to use the estimator Y^f^J = E E ^g5 +177 E ^gs- 

1=0 k=l k=l 


(4.21) 

gs-nv ■ Of course. 


the bias of this estimator is given by the bias of the Ninomiya-Victoir scheme. Thanks to its 
weak order 2, we hope to decrease the value of L, and so to reduce the computation time. We 
can also use the Ninomiya-Victoir scheme at each level and choose as follows 


= / X 


or 


and 


vO _ 
^NV — 


{f +f {x^’^’-^)) 


(4.22) 

(4.23) 


7^ — 

^NV — ~ 


X 


Ary,277? , r ( yNV,2\-p 




(4.24) 


+ f[X^ 

'xr^‘-^’^)+/(xr^'--'^)),VZe{l,...,L}. 

Actually, there is an abuse of notation in (4.24), we use the same notation rj for the 2^-dimensional 
vector (r/i,..., 772*) of the independent and identically distributed Rademacher random variables 
needed to generate the Ninomiya-Victoir scheme on the fine grid with 2^ steps and for the 
2*“^-dimensional subvector ( 7 / 1,773 ... used to generate the Ninomiya-Victoir scheme on 

the coarse grid with 2^~^ steps. The extraction of the 2^“^-dimensional vector from the 2^- 
dimensional one is aimed at reducing the variance. As previously, we obtain the same rates a 
and (3, but the main drawback is the simulation of six schemes at each level I G {1,...,L — 1} 
instead of three. Reasoning like in the proof of Proposition 4.3, since 

Z^NV = Zhs-NV+f (f (x^^’^‘-')-f (X" 


one obtains: 
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Proposition 4.4 Assume that f E with bounded first and second order deriva¬ 
tives, b £ (M"',]R"') with bounded first and second order derivatives, G G 

{1,... ,d}, with bounded first and second order derivatives and with polynomially growing third 
order derivatives, and that da^,\/j,m G {1,... ,d}, has bounded first order derivatives. Then: 


Vp > 1,3c G m;,v/ G n*, 


E 




where is defined by (4.24). 


4.2 Multilevel Richardson-Romberg extrapolation 


Recently, in [7] Lemaire and Pages developed a new method called multilevel Richardson- 
Romberg extrapolation (ML2R). This method combines the ideas behind the multilevel Monte 
Carlo approach and the multi-step Richardson-Romberg extrapolation introduced in [10]. Ac¬ 
tually, the multilevel Richardson-Romberg extrapolation can be seen as a weighted version of 
the multilevel Monte Carlo estimator. Adapting the notation of Lemaire and Pages [7], the 
multilevel Richardson-Romberg extrapolation estimator is built as follows 

= (4.25) 

1=0 k=0 

where (■^i)o<KL i<k<Mi independent random variables satisfying (4.2), (4.3) and a bias error 
expansion 


R 


3aG]R;,3RGN*,3c;,. 


,cfi G M,V/ G N,E f[Xfi 


-Y = 


(4.26) 


i=i 


where hi = T/2* is the time step. As previously, a is the order of weak convergence of the 
discretization scheme. By introducing the weights {Wi)o<i<L-: ^ smaller bias^ by 

canceling the successive bias terms in the expansion (4.26). Following [7], the computational 
complexity of Yml 2 R, denoted by Cml 2 R is defined as Cmlmc^ except that we do not take into 
account the weights {M)o<i<l- Under some assumptions (see [7] for further information), the 
optimal complexity C\^^ 2 r is given by Theorem 3.11 in [7], which states that C1^^2R depends 
on a, and the variance convergence rate ^ of denoted as previously by (3 

3c2GM+,V/gN*,v(z') <^. (4.27) 


• C*ml 2R = 0{e if /? > 1, 

• C*ml 2 R = O (e"^ log (i)) A 13 = I, 

• C*ml 2 R = O (e-2exp (^-^^2\og (2)log(i))) if ^ < 1. 


^See [7] and [10] for more details. 

^In [7], Lemaire and G. Pages assume that: VZ G N*,3Vi G R+,E 
easily adapt the proof with assumption (4.27). 



/(Tt)||' 


< Fi/if. One can 
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Similarly to the multilevel Monte Carlo method, the best complexity, obtained when /3 > 1, 
is the same as in a simple Monte Carlo method with independent and identically distributed 
unbiased random variables. With a view to achieving this complexity by applying Theorem 4.2 or 
Proposition 4.4 we will choose {^gs)o<i<l i^Nv) o<i<L Z^^fy = f . Here, we 

recall the asymptotic^ optimal parameters for the multilevel Richardson-Romberg extrapolation 
estimator: 


L* = 


- + log2 (T) 


+ - log2 
a 




(4.28) 


Mt = , 

L* 

Wi = Y.w„ 

3=1 


(4.29) 

(4.30) 


where: 


N* = [1 + 


and 


\L*-j 


Wj = (-1) 
ql<x{l + e) 




3 L*-j 

j-[ (1 _ 2-*:«) n (1 - 2"^“) 

k=l k=l 


(4.31) 




(4.32) 


E Qt = 

k 1=0 


2a {L* + 1) 


V (/ (Xt)) (^1 + 0 (^1 + f: \Wi\ (2-4^ + 2-§('-i)) V2' + 2'-i^^ 


Qo + E Qt (2' + 2'-i) 


1=1 


e = t-2 


C2 


V(/(Xt))' 


(4.33) 

(4.34) 


4.3 Numerical experiments 

In this section we present numerical tests in which we compare the multilevel Monte Carlo 
and the multilevel Richardson-Romberg estimators. Although we have not proved a theoretical 
expansion of the bias like (4.26) for the Ninomiya-Victoir and the Giles-Szpruch schemes, we 
will use these schemes in the multilevel Richardson-Romberg estimators (see [4] and [9] for 
extrapolation methods based on the Ninomiya-Victoir scheme). More precisely, we compare the 
following estimators: 


The multilevel Monte Carlo estimator with the Giles-Szpruch scheme 


yGS 

^MLMC 


7l,k 

"GS 


1=0 ‘ k=l 

where and Zj^g are respectively given by (4.12) and (4.13). 


®When e goes to 0. 
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• The multilevel Monte Carlo estimator with the Ninomiya-Victoir scheme 


M,* 




NV 


MLMC 


Y — Yz' 

MT ' 


1=0 


l,k 

NV 


I u _ 


k=l 


where Z^y and are respectively given by (4.22) or'^ (4.23) and (4.24). 

The multilevel Monte Carlo estimator with the Giles-Szpruch scheme from level 0 to level 
L* — 1, and the coupling between the Ninomiya-Victoir and the Giles-Szpruch scheme at 
the last level L* 


vGS-nv 
^ MLMC 


L*-i , .. 

Y J7^Y^gs + TTr Y 


M, — 

1=0 I k=l 


L*,k 

M* ^ “GS-W 
k=l 


where ZQg_j^y is given by (4.19). 

The multilevel Richardson-Romberg estimator with the Giles-Szpruch scheme 


L* Af;* 

- yl,k 

ML2R — 2_^ 2^ ^GS- 

1=0 * k=l 


i>GS 


The multilevel Richardson-Romberg estimator with the Ninomiya-Victoir scheme 

L* 


VNV 
^ ML2R 


yi,k 

2^ Mr 2^ NV 

1=0 ‘ k=l 


Here, Z^y is given by (4.22). 


4.3.1 Clark-Cameron SDE 


For our hrst numerical test, we consider the Clark-Cameron SDE with drift which is dehned as 
follows 

r dUt = StdW} 

\ o (4.35) 

[ dSt = fidt + dW^ 

where /z G M. In this 2—dimensional stochastic differential equation, the diffusion coefficients 


are given by 
Stratonovich drift is given by 


M and the drift coefficient is 6 f ^ 
1/ Vs 


a 


u 


= {6-i(a^V+8aV)}('“) 


0 1 
0 0 


/O 0 

0/ lo 0 


The 


These functions are smooth and satisfy the assumptions of Theorems 2.3, 3.2, Propositions 4.3 
and 4.4. By a straightforward calculation, the Giles-Szpruch scheme is given by 


UtZ = 




GS , cGS I 


ife+1 




-V 


1 


wZ - z 


K,. - Z, 


^ Cl = C++ (Cl - c) 


(4.36) 


^The choice of level 0 will be discussed later. 
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and the Ninomiya-Victoir scheme is given by 




NV,ri _ 


ik + l 


f/.r-’' + - wl) + -m(4+i - k) (wl^, - wl) 

+ hm».n (w'h. - K) (K» - K) 


(4.37) 


= s<r'’+«/■+! - </■)+- K) ■ 


Before comparing these estimators, we will illnstrate Theorems 2.3, 3.2, Propositions 4.3 and 4.4. 
In order to check the strong convergence rate of the Ninomiya-Victoir scheme, we will look at the 
expectation of the square L^-norm of the difference, at time T, between the schemes with steps hi, 
and hi-i, simulated with the same Brownian path. Denoting by 

and it follows, from Theorem 2.3 that 


E 


vNV, 2‘,V vNV,2^-^,ri 

rjn rjn 




(4.38) 


For the simulations, we choose the initial conditions C/q = Vq = 0, the final time T = 1 and the 
parameter /i = 1 . 


In figure 1, the blue line shows the behavior of log 2 ( E 


vNV,2hr] ^GS,2'-,ri 

./». rj-i y^ / 


yNV,2^,71 yNV,2‘-\7i 

yY rp y\.rp 


and the red 


line shows the behavior of log 2 ( E 


rGS,2h 


as a function of the discretization 


VIonte Carlo method with Mi = 10® 


level 1. These expectations are estimated with a standard 
samples for all 1. This choice ensures that the confidence intervals are very tight, that is why 
they are not represented in our plot. The blue line illustrates the strong convergence order 
of the Ninomiya-Victoir scheme. As expected, we obtain a line with slope -1. The red line 
illustrates the strong convergence order of the coupling between the Ninomiya-Victoir and the 
Giles-Szpruch scheme. It follows, from Theorem 3.2 that 


E 


vNV,2\ri 

Va. rp 


-X, 


GS,2‘ 


< 


221- 


(4.39) 


Again, as expected, we obtain a line with slope -2. These numerical results are consistent with 
Theorems 2.3 and 3.2 stated and proved in this paper. 

To illustrate Propositions 4.3 and 4.4, we choose a smooth payoff function, satisfying the as¬ 
sumptions of Propositions 4.3 and 4.4: f {u,s) = cos(ri). In figure 2, the top plot shows the 

i^GS-Nv)'^ ) defined by (4.19) whereas the bottom plot shows the behavior 


behavior of log 2 (^E 
of log 2 (^\{Z^NV? 


defined by (4.24). Both lines have slope -2. 

By increasing the value of fi, we noticed that the theoretical rate of convergence is reached 
for larger and larger values of 1. For small values of I, the variance decreases faster than the 
theoretical rate. Figure 3 shows this phenomenon for Zj^^y, with the payoff f{u,s) = . 

Actually, by choosing this payoff we can check that 


E 




NV 


= 2 


-4Z 




+ 2 - 3 ^ ( + 


11 


2rrb 


1545, 
512 ' 


+ 2 


-21 


163 ^4 


1024' 


(4.40) 


The details of this tedious calculation are postponed to the Appendix. The previous formula 
(4.40) contains higher order terms which overshadow the theoretical behavior of the variance. 
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Figure 1: Strong convergence order. Strong error (y-axis log 2 scale) as a function of I (x-axis) 



Figure 2: Variance convergence order with f{u,s) = cos(u). Second order moment (y-axis log 2 
scale) as a function of I (x-axis). 


The following plot shows the behavior of log 2 (E 


^ and for small values of I, the ratio E 
the leading term is 2“^^ Asymptotical 


(^w)' 


i^Nv) 


E 


i^Nv) 


as a function of 1. For large values of 
2 


is close to 16, which shows that 


y, the slopes of the curves are 2. From a numerical point 
of view and given the structure of multilevel methods, this is an important point to emphasize. 
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In particular, the choice (4.29) of parameters iii the multilevel Richardson-Romberg 

estimator is based on asymptotic properties and will not be optimal when this asymptotic 
behavior fails for the first levels. 



Figure 3: Variance convergence order with f{u,s) = . Second order moment (y-axis log 2 

scale) as a function of I (x-axis). 


Now we present the practical procedure used to implement the multilevel estimators. Putting 
together the elements already discussed, the algorithm that we use for the multilevel Monte 
Carlo with the Ninomiya-Victoir scheme or the Giles-Szpruch scheme is as follows. We begin by 
estimating the weak error constant ci in (4.7), the constant C 2 which comes from the variance 
estimation (4.8) and checking the orders of weak and strong convergence. When the asymptotic 
behavior (4.7) of the bias of the scheme is satisfied, one has 


E 



Cl (1 - 2“) 

rvj - 

2^1 


(4.41) 


Using a regression with few values of {l, |E [Z^] |), we estimate ci and check the order a of 
weak convergence. In the same way, we estimate C 2 and check the strong order /3 of variance 
convergence to 0, using a regression in (4.8). Then we estimate V (Z^) using a standard Monte 
Carlo estimator Vq. After that, for a given e we define L* using (4.9) then we set 


(4.42) 



and 


2 I 5^ 
e2V Ai2h/3+i) 



(4.43) 


When we use the Ninomiya-Victoir scheme we have the choice between Z^y = f and 

Z^v = 5 . The second choice reduces the variance of level 0 if 
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gffgg^jyg^y depends on rj. So, in general, using Z^y = 5 
reduces the sample size of the multilevel Monte Carlo estimator. Thus, although we use two 
schemes in the level 0 , the method is slightly faster with this choice in practice. As already 
mentioned, for the Giles-Szpruch scheme we choose Aq = 1 and A; = 5/2, V/ € {1, • • • ,L*}, to 
balance the lower cost of level I = 0. Following this idea, for the Ninomiya-Victoir scheme we 
choose the same sequence if Z^y = ^ , and we propose to choose 

Ao = 1 and Xi = 5,V/ G {1,... ,L*}, if Z%y = f 

Let us discuss the implementation of the multilevel Monte Carlo estimator with the Giles- 
Szpruch scheme from level 0 to level L* — 1 and the coupling between the Ninomiya-Victoir and 
the Giles-Szpruch scheme at the last level L*. The practical procedure is slightly different. As 
already discussed, in the case of given by the bias of the Ninomiya-Victoir 

scheme, so we begin with the estimation of the weak error constant ci using the Ninomiya- 
Victoir scheme. The next step is to estimate the constant C 2 using the Giles-Szpruch scheme. 
Then, we estimate V {ZQg) (respectively V {ZQ*g_]yy)) using a standard Monte Garlo estimator 
Vq^ (respectively Vgs-nv)- Finally, we define L* using (4.9) and set 


K = 


'vSs 


'^oV^g + 


L*-l 

E 

i=i 




-NV 


(4.44) 


VZg {!,...,= 


(\A^+ E 


Ml* = 


2 


L*-l 


GS-NV 


Xl*2L* 


\/M&+ /iXFiV + yVFvTTi; 


GS-NV 


i=i 


(4.45) 
. (4.46) 


We suggest to choose Aq = 1, A/ = 5/2, G {1,..., L* — 1} , and A^* = 9/2 to balance the 
higher cost of level L*. 


Since all parameters are explicit, implementing the Multilevel Richardson-Romberg estimator is 
quite simple. As noted in [7], we only need to estimate V (/ (X^)) and the constant C 2 in (4.8) 
which comes from the variance estimation. The variance V (/ (X-r)) is estimated using a crude 
Monte Garlo method. 


Now we present our numerical tests in which we compare the computing time of each estimator 
as a function of the upper bound, denoted by e, on the root mean squared error. For our first 
test we choose a smooth payoff f{u, s) = cos(tt). We estimate the two constants ci and C 2 using 
the above-mentioned procedure. To compute our regression, we estimate E [Z*] and V [Z^] for 
Z G {1,... ,4}, using a standard Monte Garlo method. The sample size used must be adjusted 
to get a rather good estimate, but without spending too much time during this step. In our 
numerical experiment, we choose a sample size M = 10^. Using this approach, we estimate 
the theoretical values of the orders of weak and variance convergences. More precisely we get 
a = 1, /3 = 2 for the Giles-Szpruch scheme and a = 2, /3 = 2 for the Ninomiya-Victoir scheme. 
In figure 4 is depicted the GPU-time in seconds (in log 2 scale) of each multilevel method as a 
function of e (in log 2 scale). It provides a direct comparison of the performance of the different 
estimators. The red line is for yMLMC^- This line is below the other lines, which indicates clearly 
that, for this experiment, is faster than the other estimators. Moreover, we observe a 

rather close behavior of the Multilevel Richardson-Romberg estimator and the Multilevel Monte 
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Carlo estimator. Indeed the black line, representing Y^Ymc close to the black dashed line 
representing Y^Y 2 R- Similarly the blue line, representing is close to the blue dashed 

line representing Y^^ 2 r- Finally, one can notice that all slopes are equal to —2, which indicates 
that all these estimators achieve a O (e~^) complexity. 



Figure 4: Clarl-Cameron SDE with f{u,s) = cos{u), CPU-time in second (y-axis log 2 scale) as 
a function of e (x-axis log 2 scale). 


To measure the efficiency of with respect to other estimators, we plot in figure 5 the 

following CPU-time ratios: 


CPU - time Y 


R = 


CPU-time 


(4.47) 


The estimator is about 1.1 to 1.6 faster than or Y^^ 2 r when e goes from 2 

to 2“’^. In comparison with Y^^a^ i '^mlmc '^ml 2 R perform poorly. 

In order to understand what is going on, let us provide a theoretical calculation of the CPU-time 
for the multilevel Monte Carlo estimator. Denoting by r; the theoretical computing time of level 
I £ {0 ,..., L*}, one has 

r' oc M;2'. (4.48) 
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0 ^ 1.6 



Figure 5: Clarl-Cameron SDE with f{u,s) = cos{u), CPU-time ratios (y-axis) as a function of 
e (x-axis log 2 scale). 


Replacing^ M*, one can write 

r' = Cl (e) 2"^(^)2' = Ci (e) 2"'(^) (4.49) 

The theoretical computing time, denoted by r, is given by 

L*(e) 

r(e) = ^ T; = ^ Q(e)2”^(^). (4.50) 

1=0 1=0 


In the multilevel Monte Carlo estimator studied in this paper Ci = Ci,\/l G {1,..., L* — 1}, 
then one has 


r(e) = 


Co (e) + ( 2 -^ - 2 -^‘«^) + Cl* (e) if/3 ^ 


1 


Co (f) + (C* (e) — 2 ) Cl (e) -|- Cl* (c) 


if /3 = 1. 


(4.51) 


Now, it is easy to understand why Y^lm^ faster than Ym^^mc- ^ matter of fact, the 
two estimators are very close, and in our numerical experiments we observe that C^^ (e) ~ 
(jGS-NV Since using a scheme with second order of weak convergence provides a lower 
optimal last level L*, in view of (4.51), we understand why, in general, we can state that 


For Ymlmc and YmYmc^ the constants 


®Obviously, the constant C (e) depend on the estimator, 
are given by formulas (4.42) and (4.43): Co (e) = ( VAoU° + X) \/c2Aj2fF“^) ] and Ci (e) = 


L* 

E 

j=i 


(4.45) and (4.46). 


V AoU° + E \/c 2 AjI, VZ € {1,.. ., I/*}. For the constans are given by formulas (4.44), 

1=1 J 
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'^ffhMC — poor performance of Ymimc ^ml 2 R reflects the use of six schemes 

in Z^j^y. 

For our second experiment, we only change the payoff. We choose the non-smooth payoff 
f{u,s) = tt+. Theorem 5.2 in [6] gives the lower bound /3 = 3/2 for the Giles-Szpruch scheme. 
Their proof is, in some ways, generic and it can easily be adapted to the Ninomiya-Victoir 
scheme. This is enough to keep the O complexity. To determine the actual values of /3 and 
a, we rely on the numerical results. Using the same automatic process, we get for the Ninomiya- 
Victoir scheme a = 3/2 and (3 = 3/2. The non-regularity of the payoff affects both the weak 
and the variance convergence rates. With regard to the Giles-Szpruch scheme, the regression 
procedure leads to a = 1 and /3 = 2, but the situation is quite confusing. Indeed, we noticed 
that the asymptotic rate /3 = 3/2 is reached for / > / = 5. Figure 6 illustrates this inflection. 
The blue line is the estimation of (V whereas the red line is the regression on the 

first four values. The two lines diverge at level I = 5, which show clearly the inflection. 



Figure 6: Clarl-Cameron SDE with f{u,s) = Variance of the Giles-Szpruch scheme (y-axis 
log 2 scale) as a function of I (x-axis). 


Here, assigning a value for (/3,C2) to implement Ymlmc^ ^ml 2 R using respec¬ 

tively (4.42)-(4.43), (4.44) to (4.46), and (4.29) to (4.34) may not be convenient. We suggest to 
apply the numerical procedure described in the following remark to implement the multilevel 
estimators. 


Remark 4.5 In the case of the Clark-Cameron SDE with = l,Uo = 0, Sq = 0 and for 
a smooth payoff, everything is going as expected, but in some cases (see the Heston model or 
the Clark-Cameron SDE with a large fi) estimating (/3,C2) may he difficult, especially when the 
theoretical rate of convergence is reached for a level I > 2 and this may affect the efficiency of 
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the multilevel methods. To get around this problem, a reasonable criterion is to compare I and 
the last level L* (e). If L* (e) < I, we decide to use the values obtained by the regression and 
use the usual formulas^ to compute for both methods. If L* (e) > I, we estimate 

V [Z*] for I G {O, using standard Monte Carlo Method. Then, we approximate, for 

I G {r+ 1, ..., L*}, V [Z*] by where is the estimation of Y Z*” and ft is the 


theoretical order of convergence of the variance. Finally, we compute Mf fori G {0, 
using (4.10). As regards the multilevel Richardson-Romberg estimator, we do not recommend its 
use in this case. 


In our second experiment, the Giles-Szpruch scheme only appears to be problematic. Indeed 
the values of L* (e) are given by: 


e 

2^ 

2-s 

2-y 

2^ 

2^ 

2“^^ 

2-14 

2^ 

2^ 

2^ 

VGS 
^ MLMC 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

vGS-NV 
^ MLMC 

3 

4 

4 

5 

6 

6 

7 

8 

8 

9 

VGS 
^ ML2R 

3 

3 

4 

4 

4 

4 

4 

5 

5 

5 


If e G 2“^^, 2“^®}, for Ymlmc^ since I is exactly equal to L{e), we are in a borderline 

situation. Nevertheless, we keep in the following figures the performance of this estimator 
for e G {2“^^, 2“^®, 2“^®|. For the multilevel Monte Carlo estimators with the Giles-Szpruch 
scheme, we apply the modified procedure of Remark 4.5 if necessary. Figure 7 compares the 
computing time of the estimator, with the previous graphical conventions. Unlike the previous 



Figure 7: Clarl-Cameron SDE with f{u,s) = u+, CPU-time in second (y-axis log 2 scale) as a 
function of e (x-axis log 2 scale). 


®(4.42)-(4.43) for Y^Imc and (4-44) to (4.46) for YmlmJ^ and (4.29) to (4.34) for Yml 2 r. 
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experiment, the two fastest estimators are Y^lmc ■ Although we lose the second 

order of weak convergence, the estimator is about 1.3 to 2 faster than Ymlmc This 

is due to the degradation of the variance convergence order (3 from 2 to 3/2 in comparison with 
a smooth payoff. Indeed, thanks to formula (4.51), one can see that, in the multilevel Monte 
Carlo methods, all things being equal, the gain in computing time due to the introduction of 
a scheme with high order of weak convergence in the last level is all the more significant that 
(3 is small. This explains why Y^Ymc performs very well. Despite using six schemes in Z^p^y, 
^MLMC goes up to 1.1 faster than (see Figure 8). In contrast, the use of a scheme with 

high order of weak convergence like Z^pjy in the multilevel Richardson-Romberg does not appear 
to counterbalance its complexity. This difference of behavior is related to the dependence of 
L* (e) on a . In the multilevel Monte Carlo estimators, the dependence is of the form 1/a (see 
(4.9)) which provides better results as alpha increases than the multilevel Richardson-Romberg 
estimator where the dependence on a is given by (4.28). 



Figure 8: Clarl-Cameron SDE with f{u,s) = u+, CPU-time ratios (y-axis) as a function of e 
(x-axis log 2 scale). 


4.3.2 Heston model 

The Heston model is an asset price model which assumes that volatility, denoted by V, evolves 
according to an autonomous Cox-Ingersoll-Ross SDE: 

f dUt = {r-^Vt)dt + ,/VtdWl 

[ dVt = k{ 9 - Vt)dt + ay/VtdWi- 

The asset price S is given by St = exp(Uf). We assume, for simplicity, no correlation between 
the Brownian motion driving the asset price and the volatility process. We also assume that 


34 





























2k6 > to ensure that the zero boundary is not attainable for the volatility process. The 
main difficulty is located in 0, where the square root is not Lipschitz. In this 2—dimensional 

( ^/u^ 9 / 0 

.v) \aJv, 


model, the diffusion coefficients are given by 


0 


o 


and the drift 


coefficient is h 


r — iu 
k{ 9 — v) 


. The Stratonovich drift is given by 




r- 


r — 2 ^ 
k{ 6 — v) 



^K{d-V)-^ 

Then, the Giles-Szpruch scheme is given by 

' vSi = vlf +«(«- Vif) '•+- * 


= f'.f + (>■ - 2 


\v,T) ft + 


(4.53) 


Setting ^ = 9 — the Ninomiya-Victoir scheme is given by 






NV,v 


If 1 A , 1 

2V 2^/ 2k 


Atf - C) (exp 


- 1 


^NV,v ^ 

*fc+§ 


TrNV,T] _ TrNV,ri 

^T+§ “^T +1 


(4.54) 


u, 


NV,rt 


tk+l 


= u, 


NV,v 


^k+l 


+ 


r - -C] h + 


2k 


V, 


NV,ri 


*fc +2 


— ^ ) ( exp ( --Kh ) — 1 


In these formulas close to our implementation of the scheme, the evolution from 

respectively , to > respectively corresponds 

to the integration of the ODE directed by the vector field on half a time step whereas 
the evolution from (to (corresponds to the integration of the 

\ t' _L T. I? _L _ / \ !> _L _ Zr*_L _ / 


fc+^ / \ k-\-4 

Brownian vector fields. The Giles-Szpruch scheme and usual schemes such as the Euler scheme 
are not well defined since they can lead to negative values of the volatility process for which 
the square root is not defined at the next step. Assuming ^ > 0, the Ninomiya-Victoir scheme 
is well dehned and the volatility process is always positive (see [1]). For .^ < 0, in section 3.1 
of [1], Alfonsi proposed a modification of the Ninomiya-Victoir scheme preserving the positivity 
of the volatility and the weak order two. For the simulation studies, we choose, as in [6], 
5'o = Do = l,r = 0.05, T = 1, k = 0.5,9 = 0.9 and a = 0.05. Then ^ = 0.89875, so that the 
Ninomiya-Victoir scheme is well defined. Using this parameters, we do not observe negative 
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values for the volatility with the Giles-Szpruch scheme. We choose to price the at the money 

call option. This corresponds to the payoff f{u,v) = exp {—rT) (exp {u) — 1 Estimating the 

multilevel parameters, we obtain, a = 1 and /3 = 2 and not 3/2 as predicted by the analysis. For 

the Nynomiya-Victoir scheme the estimation of (a, ci) leads to (2, —3.2 x 10“^). Since ci is very 

small, formula (4.9) can lead to negative values. If this occurs we set L* (e) = 1. We also observe 

that V decreases very quickly and faster than the theoretical rate for the first levels. 

Actually, an analogy can be drawn between the Ninomiya-Victoir scheme for the Clark-Cameron 

SDE and the Ninomiya-Victoir scheme for the Heston SDE since we have the same structure. As 

/ \ 2 


a matter of fact, when 




rewrites: 



This equation, similar to + // (tfc+i — tk) + ~ ™ Clark-Cameron 

SDE, is the only place where the Brownian increment AVE^ appears in the Ninomiya-Victoir 
scheme for the Heston model. In the dynamics of the U component the Brownian increment AW ^ 

is multiplied by in the Clark-Cameron SDE and 

the Heston SDE. Then, the presence of a non-zero drift probably explains the existence of the 
higher order terms that disrupt the theoretical behavior like in formula (4.40). Figure 9 illustrate 
this phenomenon. The blue line is the estimation of (V whereas the red line is the 

regression on the first four values. To be precise, we estimate V for / G {I,..., 7} using 

M = 10^ samples to get pretty good estimations, but in practice, M = 10® would be enough to 
implement the multilevel estimators. The regression leads to ,0 = 3. As in the Clark-Cameron 
SDE with f{u, s) = s_|_, the two lines diverge at level I = 5. 

So to implement the multilevel estimators, we compare I and L* (e) as already mentioned in 
Remark 4.5. The values of L* (e) are given by 


e 


2-8 

2-y 

2^ 

2^ 

2-12 

2-13 

2^ 

2^ 

2-ib 

VNV 
^ MLMC 

1 

1 

1 

1 

1 

1 
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2 
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3 

vNV 
^ ML2R 

2 

2 

2 

2 

3 

3 

3 

3 

3 
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We notice that even if e is very small, L{e) < 1. So we can implement the multilevel estimators 
using the standard automatic procedure. With regard to the Ninomiya-Victoir scheme, we 
remark that: V + f « V ^, so we naturally decide 

to implement the multilevel Monte Carlo with Z^y = / . In the following plots we 

compare the hve estimators. This time, the fastest estimator is YmYjIR- This is due to the 
outstanding value of (3 observed for the Ninomiya-Victoir scheme. The poor performance of 
'^MLMC explained by the high variance at the level 0, while the variance at the higher levels 
are very small since the numerical value of /? is 3. 

Figure 11, which represents our CPU-time ratios defined as previously, emphasizes that YJ^Y, 2 R 
is about 1.75 faster than Y^Y 2 R faster than Y^^Ymc^ since the black dashed 

curve is always below 1. 
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Figure 9: Heston SDE with f{u,v) = exp (—rT) (exp (u) — 1 Variance of the Ninomiya- 
Victoir scheme (y-axis log 2 scale) as a function of I (x-axis). 



Figure 10: Heston SDE with f{u,v) = exp(—rT) (exp (tt) — 1 CPU-time in second (y-axis 
log 2 scale) as a function of e (x-axis log 2 scale). 
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Figure 11: Heston SDE with f{u, v) = exp {—rT) (exp [u) — 1 CPU-time ratios (y-axis) as a 
function of e (x-axis log 2 scale). 


5 Conclusion 


In this paper, we have improved the multilevel Monte Carlo estimator of Giles and Szpruch [6] 
using a coupling between the Giles-Szpruch and Ninomiya-Victoir schemes at the last level of 
the MLMC estimator, which generalize their antithetic method. When the payoff is Lipschitz 
and piecewise smooth, which is very common in finance for example, the gain is amplihed since 
(3 = 3/2. We have also highlighted a strange phenomenon: sometimes the numerical rate of 
convergence of the variance can be better than the theoretical one, at least for the levels used in 
the multilevel methods. This illustrates the presence of higher order terms which overshadows 
the theoretical behavior. Therefore, we emphasize that the estimation of the rate /3 and its 
associated constant C 2 should be done cautiously, since the optimal parameters of the multilevel 
estimators depend on them. 


6 Appendix 

Let I G N* and N = 2^~^. We recall that the Glark-Gameron SDE is defined as follows 

r dUt = Stdw} 

1 dSt = fJ-dt + dW^ 
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and the Ninomiya-Victoir is given by 
if Vt^+i = -1 


' C? = t'.r-”+5^’” (<*. - K) + 5^“ (‘-=+1 - *-=) (<.« - K) 

(<« - K) {K,^ - K) («'2) 

. C;’ = f (*‘+i - *'=) + - K) 

if Vtk+i = 1 

_' Kif= t'.r-’ +•sr’” (<.*. - w'i)+5^“ (‘*+1 - ‘‘) 

. C:’ = sr^’+f (**+■ - ‘‘) + (<„ - M'i) ■ 

We define and as 1 , i respectively. 


Remark 6.1 does not depend on p, so 


The evolntion of yNV, 2 '- i,? 7 ^ coarse grid with time step T/2^ ^ is given by 


fjNV,2^-^ ,ri _ j-jNV,2‘-^,ri , nNV,2^-^ ,ri 


n + o/^ (4+1 - 4) < 


(6.4) 


+ 5 (<„ - <) (<„ - W'- 

i CW" = <’■"'■■ ’+M (‘‘«- *‘)+('^L. - o ^ 

Similarly, the evolution of yNV, 2 ’-,ri\ £j^g with time step T/2^ is given by 




+ _ Vb/ - W 






fc+^ 


(6.5) 




and: 




^*l+. - K.. 




*'=+i 


<« - K,i 


C?‘’ = CT"+" (‘‘+' - ‘‘+^) + - K.i) ■ 


( 6 . 6 ) 
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By a straightforward calculation, we get 

+ itt+i - tt) (wl^, - ir■ ^.) + If, (tf+. - ff) (n/,y. - < 




+ l *fc+i 


+ 


Wf 








n 


+ 




w/..- w, 


tk+1 




. c ?= sr '''” +'• (<<,+1 - <0 + ( w , i ^ - wi ) ^ 


(6.7) 


The antithetic scheme is defined by swapping VTij. ^ — Wt ^ and Wt ^ — 


K-- 


+ {K,i - <) ('»'■.« - K,i 


W/..,- W/ 


ife+l 


+ _ Wl _ W 










-W,l]+-( PT.i,, -Wf 


tfc+i 


*fe+-j 




w/..- ITf 


ife+l 




4Tf •” = 47’"’’+M (‘^+1 - + (»■,=„ - •r.i) 

Now we define: 


( 6 . 8 ) 


^NV,2‘,r) 1 /f7W,2',, 

tfc 2 \ ^^+1 ^ 


'jNV,2\ri 


fe + 1 


and 


^NV,2\ri _ji fgNV,2\r) , ^W,2*,r)\ 

Performing a straightforward calculation, we obtain 

' +5,7"'"'- o+- 

+ 5 (<« - <) - <vi 


K,, - K 


(6.9) 


Cf'’ = C''-^+ K (4+. - *.) + -K)- 

Then, using forward induction on k, we easily get that VA: S |0,..., — l| 


j=jNV,2\r, _ fjNV,2^-\v _ n 
BNV,2kv QNV,2^-^,ri _ „ 

‘^tk+i ‘^tk+i - U- 


( 6 . 10 ) 


We want to calculate 


Y = E 


^NV 


( 6 . 11 ) 
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where 


yi _ 
^NV — 




- X K 




+ t/. 




( 6 . 12 ) 


Using 


^ (x^ + i + u;^) = Q (x + y + M - Q (z + 


and (6.10), we get 


+ ^{^x-y-u-vf + ^{^y-x-u-vf 

+ ^{^u-x-y-vf + ^{^v-x-y-uf 

1 / \2 
--[z-w) 


1 


(6.13) 


y = —E fz^l 
16 L J 


(6.14) 


where 
Z = 




NV,2 

T 


- u. 


W,2*-l,»? _ jjNV,2‘-^-rj 


Zo = 


Z3 = 
Z4 = 


(6.15) 

(6.16) 

(6.17) 

(6.18) 

(6.19) 

( 6 . 20 ) 

In order to get an explicit expression for Zj, i G {0,..., 4}, we compute the following differences: 

<?'•” - <?'■’ = (y,,. - k) iK,i - K) 

/ X / ^ X ^ ( 6 - 21 ) 


To lighten the previous expression, we introduce 

Zi = 1 (sc"''-"'-’ - c"''-"'’” - 1 ;"''-"- c/"*'’"'’-” 

1 (st/"''-"''-’' - 

- (3U^^’‘^^~^ - _ jjNV,2\rj _ jjNV,2\-ri 

T T T T 

- (3U^^’‘^^~^ - _ jjNV,2 \-tj _ jjNV,2\’q 

4 V T T T 


and 
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-»»« {K,i - <) {K,i - '^‘ 

<r- c?'"=f'.r"'’ - +f ^ (<„ - 21^.1,. + n) 


( 6 . 22 ) 


+ I 


^fc+i 




w,",. - W, 


tk+1 




+ i - i.h) (K,i - K) (K^i - 
+ i (%.+j - m,i) (<„ - Wlj - 1^2^. 


cr- <?■ ■’=- ^r '-"+f ^ (<„ - ^K. + < 

+ (<« - M'.L J - '^‘l) - (K,i - 

+ ^ (%+i-i-=+i) (Ki-^‘>) ('^‘li-^‘ 

+ 5 (-?<=+! - %+l) ( 


(6.23) 




fc+4 


^4.. - 




(6.24) 


<r" - +f § 




- W, 


*fc+^ 


-w, 


^k 




k+4 


(6.25) 


and 


<?'■’ - C?‘’ = - ^.7’"'”' + f ^ {k,. - + w'i 




- W, 


*fc+^ 


- W, 


^k 




Hence, we obtain by summation 


Af-l 


4 


where 


z. = E 

k=0 

4 = (iri^ - wl) {w,l^^ - wl) 


^+1. 

(6.26) 


(6.27) 

(6.28) 
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and 


4 = - (i - i - %«)) - wQ (wl, - w. 


*fc + l *i:+l 


Then, one can express Z'^ as 




- (<« - 2^.4 + K) + («'h4 - '»'<.) ("'^4 - 

+ {k» - K^) {k» - K^i) ■ 


^N-1 


^N-l 


/N-l 


/N-1 


/N-l 


= 




.k=0 


> fe =0 


■.k=0 


■.k=0 


. fc =0 


(6.29) 


(6.30) 


(6.31) 


(6.32) 


i X] ^ X] “*"^ 2 ^ + ^14 + 44 “ 44 ^ 

V fc=o (fe,0eAjv / 


(6.33) 


where A^r = |(A:, /) £ {0,... , A — 1}^ , A: < /|. 
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Preliminary calculus: 


We begin by writing in generic form: 

4 = 4 (ir,L. - + A [w^ - w.\) 

+ 4 - 2 ir‘+ ir') + 4 (w^ - ir-) 


(6.34) 


where 



First, let us look at the expectation of zlzi: 


E 


r • '1 

r2 r 

^k^k 

=-^E 

4A^2 [ 


+ pwi + 


T 


:E 


iklk 


(6.35) 


Then, for i = j = 0 


For i = j e {1,... ,4} 


For j = 0 and z G {1,..., 4} 


For (z, j) G {(1,2), (3,4)} 


For (i,j) G 1(1,3), (2,4)1 


For (z, j) G {(1,4), (2,3)1 


E 



rj^2 

W' 


E 



IIT^ 

64iV2 + 16JV3- 


E [zlzt] = 0. 


E 


44 


7T2 

64A^2 “ i6Ar3- 



5 r 2 ^2^-3 

64A^2 + i6jv3- 


E 


44 


92^2 ^ 22^3 

64A^2 “ i6iv3' 


(6.36) 

(6.37) 

(6.38) 

(6.39) 

(6.40) 

(6.41) 


44 
































Now we look at the expectation of [z^) ■ 

92^4 


E 


( 4 )‘ 


+ 


16A^4 

rpA 


E 


(4) 


(4)^ 


+ ( 4 ) ( 4 ) + i ^ k ) (^; 


16A^4 


E 


+ 4E 


3T^ f 

+ ^^1 E 


16A^4 


( 4 )^ + UY {PlY + ( 4 )^ ( 4 )' + ( 4 )' HY 


(«fe4 + 4^1) (44 + 44 ) + ( 









+ 4E 


44 + f^kPk) i^k^i + ^k^i 



hY ( (4) + (4) + (4) + (4) - 444 - 444 


+ (4) + (4) + (4) “ 40:^5^ - 


ihi 44 + 44 + 44 + ^14 - 244 - 244 - 244 - 




+ ( 4 ) + ( 4 ) + ( 4 ) ) 
ihl (44 + 44 + 44 + ^fe 4 


nrp2 

+ ^E 


HY hi , 


Then, by straightforward calculation, for i = j = 0 


E 


0^4 


( 4 ) 


92^4 


For i = j S {1,... ,4} 


E 


(4)^ 


23ir^ 33/i2r5 34^® 


1024Ar4 256A^5 256A^6 ' 


For 2=0 and i € {1,..., 4} 


E 


(4)44)' 


37r4 , 

64iV4 + i6jv5' 


For (i, 2 ) € {(1,2), (3,4)} 


E 


(4)' 


i27r4 254 t 5 34r6 


1024iV4 256iV5 256iV6' 



(6.42) 


(6.43) 


(6.44) 


(6.45) 


(6.46) 
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For {i,j) G {(1,3), (2,4)} 


(4) 


2 


21//2r5 


E 


J ~3 
^k^k 


1024iV4 256iV5 256iV6' 


159r^ 29/r2T5 ^ 


E 

For (i,j) G 1(1,4), (2,3)1 

Now, we define: 
and 

bki = zizi + 2:^4 + 44 + 44 “ 44- 

Then ^2 can be expressed as 
N-l 

Ofc + 2 a^ai + 4 + 4 


1024iV4 256iV5 256iV6- 

a;. = (4)' + (4)%(4)%(4)'-(4f 


cij bki. 


k=0 

Observing that 


N 


E [ajbki] = 0 


(fc,z)eA^je[0;Af-iI 


we obtain 


E 


N-l 



= E 

4i 



(/c,/)GAiv 


E[z 2] = ^E[ai] +2 ^ E[afcaz] + 4 ^ E[b 

k=0 (fc,0eAiv (fc,06Ajv 


For k ^ I, by independence 


E [akai] = E [ofc] E [ai] = (E [a^])^ 


(6.47) 

(6.48) 

(6.49) 

(6.50) 

(6.51) 

(6.52) 

(6.53) 

(6.54) 

(6.55) 


this leads to 


E [Z^] = iVE [al] +N{N-1){E [ak]f + 2N{N-1)E [bh] . 
To achieve our goal, we calculate the following expectations: 


it 2 


E [ofc] = 4E (4) 


-E 


0^2 


( 4 ) 


5r2 


41V3 16Ai2 


(6.56) 

(6.57) 


E[4] 


4E 


(4)^ 


+ E 


(4)^ 


(4)' (4)^ 


+ 4 E 


{ 4 )\ 4 y 


+ E 


{ 4 )\ 4 y 


+ E 


{ 4 )\ 4 y 


, 19/r2T5 ^ 427r^ 

16iV6 + i6iv5 + 64A^4 

(6.58) 


E [4z] 



/i^T^ , /i22.5 ^ 092^4 

16A^6 + 4jv5 + 25^' 


8 (E [4^0])' + 4 ((E [zlzl]f + (E [44])' + (E [ 44 ])') 


(6.59) 
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Combining our results, we obtain 


E [Z^] = 


1 


) + 


2t^5 


V 16 

Then replacing N by 2*“^ 


16 


1 /II 


iV3 \ 32 


1 + 


1545, 


256 


1 /163, 


iV2 \ 256 


E [Z^] = 2-^' (3/r® + + 2-3' + 


32 J V 64 


Dividing by 16, we get the desired result 


F = 2 


-41 


4/r» + ) + 2-=“ 


11 2 . 1545 4 , 

-,V + —T^I+2 


-21 


163^4 


1024' 


(6.60) 


(6.61) 


(6.62) 
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